## Daniel J. Mallinson
## Reproduction Script for: 
## "Growth and Gaps: A Meta-Review of Policy Diffusion Studies in the American States"
## In Policy & Politics
## 2021-02-23

rm(list = ls(all = TRUE))

#Read necessary libraries
install.packages(c("foreign", "ggplot2", "Hmisc", "plyr", "meta","metafor"))

library(foreign)
library(ggplot2)
library(Hmisc) # Needed for minor tick marks
library(plyr)
library(meta)
library(metafor)

#Create color psublue
psublue <- rgb(red=0.0, green=0.14, blue=0.4, alpha=1)


##################################################################
####################### Systematic review ########################
##################################################################

#Read data
load.studies <- read.csv("PP Study Details.csv")
studies <- load.studies
colnames(studies)

#subset on included studies
include <- studies[which(studies$exclude==FALSE),]

##################################################################
####################### Summary of inclusion/exclusion ###########
##################################################################

#Total number of models: 538
nrow(studies)

#Total number of included models: 507
nrow(include)

#Total number of articles: 185
unique <- as.data.frame(unique(include$title))
names(unique) <- "title"
unique$study <- seq(1,nrow(unique),1)

include <- merge(include, unique, by="title", all.x=TRUE)

unique <- include[order(include$study, decreasing=TRUE),]
unique <- unique[!duplicated(unique$study),]
nrow(unique) #185

#######################################################################
################ Visualize major topics - Figure 2 ####################
#######################################################################

#Only included studies

major.topic <- table(include$majortopic)
mt.prop <- prop.table(major.topic)
mt.prop

a <- matrix(table(include$majortopic))
a <- rbind(a, 0, 0, 0,0,0)
names <- rbind(1,2,3,4,5,6,7,8,10,12,13,14,15,20,21,24,16,17,18,19)
a <- cbind(names,a)
a <- as.data.frame(a)
colnames(a) <- c("code", "count")
a <- a[order(a$code),]

jpeg("figure2.jpg", width=11, height=8.5,units="in", res=1000)
par(las=1)
par(mar=c(5,12,1,2))
names <- c("Fiscal (1)","Civil Rights (2)","Health (3)","Agriculture (4)","Labor (5)","Education (6)","Environment (7)","Energy (8)","Transportation (10)","Law (12)","Social Welfare (13)","Housing (14)","Commerce (15)","Defense (16)","Technology (17)","Trade (18)","International Affairs (19)" ,"Administrative (20)","Public Resource Mgmt (21)","Local Government (24)")
#bar <- barplot(a$count, horiz=TRUE, xlab="Total Number of Policies in Final Dataset", ylab="", yaxt='n', col=rainbow(16))
bar <- barplot(a$count, horiz=TRUE, ylab="", yaxt='n', col=psublue, font=2)
minor.tick(nx=5, ny=1, tick.ratio=.5)
title(xlab="Count of Models", font=2, font.sub=2)
axis(2, at=bar, labels=names, las=2, font=2)
mtext(c("n=507"), side=1, adj=0, line=3)
dev.off()

#######################################################################
############## Visualize Model Time Frames - Figure 3 #################
#######################################################################

times <- as.matrix(cbind(include$first_adopt, include$last_adopt))
#times <- times[-53:-55,]
times <- times[order(times[,1]),]
times <- cbind(times, c(1:nrow(times)))
times <- times[-(505:507),]
times <- times[-(1:3),]

modelyears <- 

jpeg("figure3.jpg", width=11, height=8.5,units="in", res=1000)
plot.new()
par(mfrow=c(1,2))
hist(modelyears, breaks=119, col=psublue, border=psublue, ylim=c(0,350), xlim=c(1900,2020), main="(a)", ylab="Count of Years in Models", xlab="Observed Year", axes=FALSE)
axis(1, at=seq(1900,2020,10), labels=seq(1900,2020,10))
axis(2, at=seq(0,350,50), labels=seq(0,350,50), las=2)
abline(v=c(seq(1900,2020,10)), lwd=1.5)
mtext(c("n=501"), side=1, adj=0, line=3)

plot.new()
plot.window(xlim=c(1900,2020), ylim=c(0,nrow(times)))
for(i in (1:nrow(times))){
segments(times[i,1],i, times[i,2], i, col=psublue, lwd=2)
}
axis(1, at=seq(1900,2020,10), labels=seq(1900,2020,10), font=1)
title(xlab="Year Included in Model", main="(b)")
title(ylab="Models Ordered by First Adoption Year", line=0)
abline(v=1960, lty=2, lwd=2, col="red")
mtext(c("n=501"), side=1, adj=0, line=3)
dev.off()

##################################################################
######################### Meta-Analysis ##########################
##################################################################

variables <- read.csv("PP Variables.csv")

#merge variables with study details

meta <- merge(variables, studies, by="model_id")

#Drop excluded studies
meta.include <- meta[which(meta$exclude==0),]

## Count of included models: 507

length(as.numeric(unique(meta.include$model_id)))

##################################################################
######################### Table 2 ################################
##################################################################
## Proportion of each variable type included per model

model <- as.factor(meta.include$model_id)

# Count of models with some ideology measure: 305
count.ideo <- ddply(meta.include, .(model), summarize, freq=sum(ideology))
nrow(subset(count.ideo, freq>0))
# Percentage of models with some ideology measure: 60 percent
nrow(subset(count.ideo, freq>0))/length(as.numeric(unique(meta.include$model_id)))

# Count of models with legislative professionalism: 126
count.legprof <- ddply(meta.include, .(model), summarize, freq=sum(profession))
nrow(subset(count.legprof, freq>0))
# Percentage of models with legislative professionalism: 25
nrow(subset(count.legprof, freq>0))/length(as.numeric(unique(meta.include$model_id)))

# Count of models with at least one slack resource measure: 350
count.slack <- ddply(meta.include, .(model), summarize, freq=sum(slack))
nrow(subset(count.slack, freq>0))
# Percentage of models with at least one slack resource measure: 69 percent
nrow(subset(count.slack, freq>0))/length(as.numeric(unique(meta.include$model_id)))

# Count of models with some political measure: 418
count.pol <- ddply(meta.include, .(model), summarize, freq=sum(political))
nrow(subset(count.pol, freq>0))
# Percentage of models with some political measure: 82 percent
nrow(subset(count.pol, freq>0))/length(as.numeric(unique(meta.include$model_id)))

# Count of studies with contextual measures: 461
count.context <- ddply(meta.include, .(model), summarize, freq=sum(context))
nrow(subset(count.context, freq>0))
# Percentage of studies with contextual measures: 91 percent
nrow(subset(count.context, freq>0))/length(as.numeric(unique(meta.include$model_id)))

# Count of studies with other regional measures: 139
count.geog <- ddply(meta.include, .(model), summarize, freq=sum(geog))
nrow(subset(count.geog, freq>0))
# Percentage of studies with other regional measures: 27 percent
nrow(subset(count.geog, freq>0))/length(as.numeric(unique(meta.include$model_id)))

#Count of studies with federal influence measure: 25
count.federal <- ddply(meta.include, .(model), summarize, freq=sum(fed))
nrow(subset(count.federal, freq>0))
# Percentage of studies with federal measures: 5 percent
nrow(subset(count.federal, freq>0))/length(as.numeric(unique(meta.include$model_id)))

#Count of studies with at least one demographic measure: 383
count.demo <- ddply(meta.include, .(model), summarize, freq=sum(demog))
nrow(subset(count.demo, freq>0))
# Percentage of studies with federal measures: 76 percent
nrow(subset(count.demo, freq>0))/length(as.numeric(unique(meta.include$model_id)))

#Count of studies with a neighbor measure: 370
count.neigh <- ddply(meta.include, .(model), summarize, freq=sum(neigh))
nrow(subset(count.neigh, freq>0))
#Percentage of studies with neighbor measure: 73 percent
nrow(subset(count.neigh, freq>0))/length(as.numeric(unique(meta.include$model_id)))

##################################################################
###################### Meta-Analysis #############################
##################################################################

##################################################################
########################## Table 3 ###############################
##################################################################

#################### Neighbor ####################################

# Create subset of neighbor data

neigh.data <- meta.include[which(meta.include$neigh==1),]
nrow(neigh.data) #401

#correct three effects types

# Drop 119 and 120 due to non-reporting of effect

neigh.data <- subset(neigh.data, neigh.data$model_id != 119)
neigh.data <- subset(neigh.data, neigh.data$model_id != 120)

neigh.data$effects_type #check for NAs

#i <- which(neigh.data$model_id==18) #not sure what these were for
#j <- "effects_type"

#Drop coefficients that are other effects types (effects_type==3)

neigh.data <- subset(neigh.data, neigh.data$effects_type < 3)
nrow(neigh.data) #389

#Check for 3s
neigh.data$effects_type

# Recode "significance" to be -1, 0, or 1
neigh.data$sig_recode <- 0
neigh.data$sig_recode
neigh.data$sig_recode[neigh.data$coefficient > 0 & neigh.data$significance > 1] <- 2
neigh.data$sig_recode
neigh.data$sig_recode[neigh.data$coefficient > 0 & neigh.data$significance <= 1] <- 1
neigh.data$sig_recode
neigh.data$sig_recode[neigh.data$coefficient < 0 & neigh.data$significance > 1] <- -2
neigh.data$sig_recode
neigh.data$sig_recode[neigh.data$coefficient < 0 & neigh.data$significance <= 1] <- -1
neigh.data$sig_recode

table(neigh.data$sig_recode)
prop.table(table(neigh.data$sig_recode))

#################### Ideological Distance #############################

gncp.data <- meta.include[which(meta.include$ideological_dist==1),]
nrow(gncp.data)  # 44

# Check for NAs
gncp.data$significance

gncp.data$sig_recode <- 0
gncp.data$sig_recode
gncp.data$sig_recode[gncp.data$coefficient > 0 & gncp.data$significance > 1] <- 2
gncp.data$sig_recode
gncp.data$sig_recode[gncp.data$coefficient > 0 & gncp.data$significance <= 1] <- 1
gncp.data$sig_recode
gncp.data$sig_recode[gncp.data$coefficient < 0 & gncp.data$significance > 1] <- -2
gncp.data$sig_recode
gncp.data$sig_recode[gncp.data$coefficient < 0 & gncp.data$significance <= 1] <- -1
gncp.data$sig_recode

table(gncp.data$sig_recode)
prop.table(table(gncp.data$sig_recode))

#################### Federal Intervention #############################

fed.data <- meta.include[which(meta.include$fed==1),]
nrow(fed.data)  # 38

# Check for NAs
fed.data$significance

fed.data$sig_recode <- 0
fed.data$sig_recode
fed.data$sig_recode[fed.data$coefficient > 0 & fed.data$significance > 1] <- 2
fed.data$sig_recode
fed.data$sig_recode[fed.data$coefficient > 0 & fed.data$significance <= 1] <- 1
fed.data$sig_recode
fed.data$sig_recode[fed.data$coefficient < 0 & fed.data$significance > 1] <- -2
fed.data$sig_recode
fed.data$sig_recode[fed.data$coefficient < 0 & fed.data$significance <= 1] <- -1
fed.data$sig_recode

table(fed.data$sig_recode)
prop.table(table(fed.data$sig_recode))

#################### Per Capita Income #############################

percap.data <- meta.include[which(meta.include$percap==1),]
nrow(percap.data)  # 196

# Check for NAs
percap.data$significance #3 NAs
percap.data <- percap.data[!is.na(percap.data$significance),] #remove legitimate NAs

percap.data$sig_recode <- 0
percap.data$sig_recode
percap.data$sig_recode[percap.data$coefficient > 0 & percap.data$significance > 1] <- 2
percap.data$sig_recode
percap.data$sig_recode[percap.data$coefficient > 0 & percap.data$significance <= 1] <- 1
percap.data$sig_recode
percap.data$sig_recode[percap.data$coefficient < 0 & percap.data$significance > 1] <- -2
percap.data$sig_recode
percap.data$sig_recode[percap.data$coefficient < 0 & percap.data$significance <= 1] <- -1
percap.data$sig_recode

table(percap.data$sig_recode)
prop.table(table(percap.data$sig_recode))

#################### Legislative Professionalism #########################

legprof.data <- meta.include[which(meta.include$profession==1),]
nrow(legprof.data)  # 131

# Check for NAs
legprof.data$significance

# Recode "significance" to be -2 (significant negative), -1 (ns negative), 0, 1 (ns positive), or 2 (significant positive)
legprof.data$sig_recode <- 0
legprof.data$sig_recode
legprof.data$sig_recode[legprof.data$coefficient > 0 & legprof.data$significance > 1] <- 2
legprof.data$sig_recode
legprof.data$sig_recode[legprof.data$coefficient > 0 & legprof.data$significance <= 1] <- 1
legprof.data$sig_recode
legprof.data$sig_recode[legprof.data$coefficient < 0 & legprof.data$significance > 1] <- -2
legprof.data$sig_recode
legprof.data$sig_recode[legprof.data$coefficient < 0 & legprof.data$significance <= 1] <- -1
legprof.data$sig_recode

table(legprof.data$sig_recode)
prop.table(table(legprof.data$sig_recode))

#################### Liberalism ####################################

# Citizen Liberalism
cit.lib <- meta.include[which(meta.include$berry_cit==1),]
nrow(cit.lib) #160

# Government Liberalism
gov.lib <- meta.include[which(meta.include$berry_gov==1),]
nrow(gov.lib) #91

## Identify coefficients that are significant (p < 0.05)

# Check for NAs
cit.lib$significance
gov.lib$significance

# Recode "significance" to be -1, 0, or 1
cit.lib$sig_recode <- 0
cit.lib$sig_recode
cit.lib$sig_recode[cit.lib$coefficient > 0 & cit.lib$significance > 1] <- 2
cit.lib$sig_recode
cit.lib$sig_recode[cit.lib$coefficient > 0 & cit.lib$significance <= 1] <- 1
cit.lib$sig_recode
cit.lib$sig_recode[cit.lib$coefficient < 0 & cit.lib$significance > 1] <- -2
cit.lib$sig_recode
cit.lib$sig_recode[cit.lib$coefficient < 0 & cit.lib$significance <= 1] <- -1
cit.lib$sig_recode

gov.lib$sig_recode <- 0
gov.lib$sig_recode
gov.lib$sig_recode[gov.lib$coefficient > 0 & gov.lib$significance > 1] <- 2
gov.lib$sig_recode
gov.lib$sig_recode[gov.lib$coefficient > 0 & gov.lib$significance <= 1] <- 1
gov.lib$sig_recode
gov.lib$sig_recode[gov.lib$coefficient < 0 & gov.lib$significance > 1] <- -2
gov.lib$sig_recode
gov.lib$sig_recode[gov.lib$coefficient < 0 & gov.lib$significance <= 1] <- -1
gov.lib$sig_recode

#Vote Count
table(cit.lib$sig_recode)
prop.table(table(cit.lib$sig_recode))

table(gov.lib$sig_recode)
prop.table(table(gov.lib$sig_recode))

#################### Initiative Available #########################

init.data <- meta.include[which(meta.include$initiative_avail==1),]
nrow(init.data)  # 33

# Check for NAs
init.data$significance

# Recode "significance" to be -1, 0, or 1
init.data$sig_recode <- 0
init.data$sig_recode
init.data$sig_recode[init.data$coefficient > 0 & init.data$significance > 1] <- 2
init.data$sig_recode
init.data$sig_recode[init.data$coefficient > 0 & init.data$significance <= 1] <- 1
init.data$sig_recode
init.data$sig_recode[init.data$coefficient < 0 & init.data$significance > 1] <- -2
init.data$sig_recode
init.data$sig_recode[init.data$coefficient < 0 & init.data$significance <= 1] <- -1
init.data$sig_recode

table(init.data$sig_recode)
prop.table(table(init.data$sig_recode))

#################### Democratic Control #########################

dem.data <- meta.include[which(meta.include$dem_control==1),]
nrow(dem.data)  # 59

# Check for NAs
dem.data$significance # 2 NAs
dem.data <- dem.data[!is.na(dem.data$significance),] #remove legitimate NAs

# Recode "significance" to be -1, 0, or 1
dem.data$sig_recode <- 0
dem.data$sig_recode
dem.data$sig_recode[dem.data$coefficient > 0 & dem.data$significance > 1] <- 2
dem.data$sig_recode
dem.data$sig_recode[dem.data$coefficient > 0 & dem.data$significance <= 1] <- 1
dem.data$sig_recode
dem.data$sig_recode[dem.data$coefficient < 0 & dem.data$significance > 1] <- -2
dem.data$sig_recode
dem.data$sig_recode[dem.data$coefficient < 0 & dem.data$significance <= 1] <- -1
dem.data$sig_recode

table(dem.data$sig_recode)
prop.table(table(dem.data$sig_recode))

#################### Republican Control #########################

rep.data <- meta.include[which(meta.include$rep_control==1),]
nrow(rep.data)  # 95

# Check for NAs
rep.data$significance 

# Recode "significance" to be -1, 0, or 1
rep.data$sig_recode <- 0
rep.data$sig_recode
rep.data$sig_recode[rep.data$coefficient > 0 & rep.data$significance > 1] <- 2
rep.data$sig_recode
rep.data$sig_recode[rep.data$coefficient > 0 & rep.data$significance <= 1] <- 1
rep.data$sig_recode
rep.data$sig_recode[rep.data$coefficient < 0 & rep.data$significance > 1] <- -2
rep.data$sig_recode
rep.data$sig_recode[rep.data$coefficient < 0 & rep.data$significance <= 1] <- -1
rep.data$sig_recode

table(rep.data$sig_recode)
prop.table(table(rep.data$sig_recode))

#################### Divided Government #########################

div.data <- meta.include[which(meta.include$divided==1),]
nrow(div.data)  # 66

# Check for NAs
div.data$significance

# Recode "significance" to be -1, 0, or 1
div.data$sig_recode <- 0
div.data$sig_recode
div.data$sig_recode[div.data$coefficient > 0 & div.data$significance > 1] <- 2
div.data$sig_recode
div.data$sig_recode[div.data$coefficient > 0 & div.data$significance <= 1] <- 1
div.data$sig_recode
div.data$sig_recode[div.data$coefficient < 0 & div.data$significance > 1] <- -2
div.data$sig_recode
div.data$sig_recode[div.data$coefficient < 0 & div.data$significance <= 1] <- -1
div.data$sig_recode

table(div.data$sig_recode)
prop.table(table(div.data$sig_recode))

#################### Slack Resources #############################

slack.data <- meta.include[which(meta.include$slack==1),]
nrow(slack.data)  # 532

# Check for NAs
slack.data$significance

# Remove 3 legitimate NAs
slack.data <- subset(slack.data, !is.na(slack.data$significance))
nrow(slack.data) #529

slack.data$sig_recode <- 0
slack.data$sig_recode
slack.data$sig_recode[slack.data$coefficient > 0 & slack.data$significance > 1] <- 1
slack.data$sig_recode
slack.data$sig_recode[slack.data$coefficient < 0 & slack.data$significance > 1] <- -1
slack.data$sig_recode

table(slack.data$sig_recode)
prop.table(table(slack.data$sig_recode))


#################### Contextual Variables #########################

context.data <- meta.include[which(meta.include$context==1),]
nrow(context.data)  # 406

## Identify coefficients that are significant (p < 0.05)

# Check for NAs
context.data$significance

# Replace NA in "significance" for studies with a t-stat greater than 1.96 (p < 0.05).
context.data$tstat[is.na(context.data$tstat)] <- .002 #Assign small value so next step works
context.data$significance[(is.na(context.data$significance)) & (abs(context.data$tstat) > 1.96)] <- 2
context.data$significance
context.data$significance[(is.na(context.data$significance)) & (abs(context.data$tstat) < 1.96)] <- 0
context.data$significance

context.data$sig_recode <- 0
context.data$sig_recode
context.data$sig_recode[context.data$coefficient > 0 & context.data$significance > 1] <- 1
context.data$sig_recode
context.data$sig_recode[context.data$coefficient < 0 & context.data$significance > 1] <- -1
context.data$sig_recode

table(context.data$sig_recode)
prop.table(table(context.data$sig_recode))

## Database Table 2
length(unique(meta$model_id[meta$neigh==1]))
length(unique(meta$model_id[meta$geog==1]))
length(unique(meta$model_id[meta$south==1]))
length(unique(meta$model_id[meta$fed==1]))
length(unique(meta$model_id[meta$ideology==1]))
length(unique(meta$model_id[meta$berry_govt==1]))
length(unique(meta$model_id[meta$berry_cit==1]))
length(unique(meta$model_id[meta$ideological_dist==1]))
length(unique(meta$model_id[meta$duration==1]))
length(unique(meta$model_id[meta$slack==1]))
length(unique(meta$model_id[meta$percap==1]))
length(unique(meta$model_id[meta$demog==1]))
length(unique(meta$model_id[meta$urban==1]))
length(unique(meta$model_id[meta$population==1]))
length(unique(meta$model_id[meta$political==1]))
length(unique(meta$model_id[meta$interest_groups==1]))
length(unique(meta$model_id[meta$initiative_avail==1]))
length(unique(meta$model_id[meta$dem_control==1]))
length(unique(meta$model_id[meta$rep_control==1]))
length(unique(meta$model_id[meta$divided==1]))
length(unique(meta$model_id[meta$profession==1]))
length(unique(meta$model_id[meta$context==1]))
length(unique(meta$model_id[meta$other==1]))
length(unique(meta$model_id[meta$intercept==1]))


##################################################################
########################## Table 4 ###############################
##################################################################

#Resource:
#http://daniellakens.blogspot.nl/2014/08/on-reproducibility-of-meta-analyses.html

## Neighbor Adoptions
# Convert odds ratios into log-odds for meta-analysis

#Break the dataset
neigh.logodds <- neigh.data[which(neigh.data$effects_type==1),]
neigh.oddsrat <- neigh.data[which(neigh.data$effects_type==2),]

#log the odds ratios
neigh.oddsrat$coefficient <- log(neigh.oddsrat$coefficient)

#Reform dataset
neigh.data <- rbind(neigh.logodds, neigh.oddsrat)

neigh.data <- neigh.data[which(neigh.data$se > 0),]
neigh.data <- neigh.data[which(neigh.data$se<100),]

# Meta-analysis
m.neigh <- metagen(coefficient, se, data=neigh.data, studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(m.neigh)

## Ideological Distance
#Break the dataset
gncp.logodds <- gncp.data[which(gncp.data$effects_type==1),]
gncp.oddsrat <- gncp.data[which(gncp.data$effects_type==2),]

#log the odds ratios
gncp.oddsrat$coefficient <- log(gncp.oddsrat$coefficient)

#Reform dataset
gncp.data <- rbind(gncp.logodds, gncp.oddsrat)

gncp.data <- gncp.data[which(gncp.data$se > 0),]
#gncp.data <- gncp.data[which(gncp.data$se<100),]

# Meta-analysis
m.gncp <- metagen(coefficient, se, data=gncp.data, studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(m.gncp)

##Federal Intervention
#Break the dataset
fed.logodds <- fed.data[which(fed.data$effects_type==1),]
fed.oddsrat <- fed.data[which(fed.data$effects_type==2),]

#log the odds ratios
fed.oddsrat$coefficient <- log(fed.oddsrat$coefficient)

#Reform dataset
fed.data <- rbind(fed.logodds, fed.oddsrat)

fed.data <- fed.data[which(fed.data$se > 0),]
#gncp.data <- gncp.data[which(gncp.data$se<100),]

# Meta-analysis
m.fed<- metagen(coefficient, se, data=fed.data, studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(m.fed)

## Per Capita Income
#Break the dataset
percap.logodds <- percap.data[which(percap.data$effects_type==1),]
percap.oddsrat <- percap.data[which(percap.data$effects_type==2),]

#log the odds ratios
percap.oddsrat$coefficient <- log(percap.oddsrat$coefficient)

#Reform dataset
percap.data <- rbind(percap.logodds, percap.oddsrat)

percap.data <- percap.data[which(percap.data$se > 0),]
#percap.data <- percap.data[which(percap.data$se<100),]

# Meta-analysis
m.percap<- metagen(coefficient, se, data=percap.data, studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(m.percap)

## Legislative Professionalism
#Break the dataset
legprof.logodds <- legprof.data[which(legprof.data$effects_type==1),]
legprof.oddsrat <- legprof.data[which(legprof.data$effects_type==2),]

#log the odds ratios
legprof.oddsrat$coefficient <- log(legprof.oddsrat$coefficient)

#Reform dataset
legprof.data <- rbind(legprof.logodds, legprof.oddsrat)

legprof.data <- legprof.data[which(legprof.data$se > 0),]
#legprof.data <- legprof.data[which(legprof.data$se<100),]

## Meta-analysis
m.legprof<- metagen(coefficient, se, data=legprof.data, studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(m.legprof)

## Liberalism
#Break the dataset
cit.lib.logodds <- cit.lib[which(cit.lib$effects_type==1),]
cit.lib.oddsrat <- cit.lib[which(cit.lib$effects_type==2),]

#log the odds ratios
cit.lib.oddsrat$coefficient <- log(cit.lib.oddsrat$coefficient)

#Reform dataset
cit.lib <- rbind(cit.lib.logodds, cit.lib.oddsrat)

cit.lib <- cit.lib[which(cit.lib$se > 0),]
#cit.lib <- cit.lib[which(cit.lib$se<100),]

## Meta-analysis
m.cit.lib <- metagen(coefficient, se, data=cit.lib, studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(m.cit.lib)

#Break the dataset
gov.lib.logodds <- gov.lib[which(gov.lib$effects_type==1),]
gov.lib.oddsrat <- gov.lib[which(gov.lib$effects_type==2),]

#log the odds ratios
gov.lib.oddsrat$coefficient <- log(gov.lib.oddsrat$coefficient)

#Reform dataset
gov.lib <- rbind(gov.lib.logodds, gov.lib.oddsrat)

gov.lib <- gov.lib[which(gov.lib$se > 0),]
#gov.lib <- gov.lib[which(gov.lib$se<100),]

# Meta-analysis
m.gov.lib<- metagen(coefficient, se, data=gov.lib, studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(m.gov.lib)

## Initiative Available
#Break the dataset
init.logodds <- init.data[which(init.data$effects_type==1),]
init.oddsrat <- init.data[which(init.data$effects_type==2),]

#log the odds ratios
init.oddsrat$coefficient <- log(init.oddsrat$coefficient)

#Reform dataset
init.data <- rbind(init.logodds, init.oddsrat)

init.data <- init.data[which(init.data$se > 0),]
#init.data <- init.data[which(init.data$se<100),]

# Meta-analysis
m.init<- metagen(coefficient, se, data=init.data, studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(m.init)

## Democratic Control of Government
#Break the dataset
dem.logodds <- dem.data[which(dem.data$effects_type==1),]
dem.oddsrat <- dem.data[which(dem.data$effects_type==2),]

#log the odds ratios
dem.oddsrat$coefficient <- log(dem.oddsrat$coefficient)

#Reform dataset
dem.data <- rbind(dem.logodds, dem.oddsrat)

dem.data <- dem.data[which(dem.data$se > 0),]
#dem.data <- dem.data[which(dem.data$se<100),]

# Meta-analysis
m.dem<- metagen(coefficient, se, data=dem.data, studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(m.dem)

## Republican Control of Government
#Break the dataset
rep.logodds <- rep.data[which(rep.data$effects_type==1),]
rep.oddsrat <- rep.data[which(rep.data$effects_type==2),]

#log the odds ratios
rep.oddsrat$coefficient <- log(rep.oddsrat$coefficient)

#Reform dataset
rep.data <- rbind(rep.logodds, rep.oddsrat)

rep.data <- rep.data[which(rep.data$se > 0),]
#rep.data <- rep.data[which(rep.data$se<100),]

# Meta-analysis
m.rep<- metagen(coefficient, se, data=rep.data, studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(m.rep)

## Divided Government
#Break the dataset
div.logodds <- div.data[which(div.data$effects_type==1),]
div.oddsrat <- div.data[which(div.data$effects_type==2),]

#log the odds ratios
div.oddsrat$coefficient <- log(div.oddsrat$coefficient)

#Reform dataset
div.data <- rbind(div.logodds, div.oddsrat)

div.data <- div.data[which(div.data$se > 0),]
#div.data <- div.data[which(div.data$se<100),]

# Meta-analysis
m.div <- metagen(coefficient, se, data=div.data, studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(m.div)

##################################################################
########################### Table 5 ##############################
##################################################################

## Neighbor Adoptions

# Regulatory policy
r.neigh <- metagen(coefficient, se, data=neigh.data[which(neigh.data$boushey==1),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(r.neigh)

mo.neigh <- metagen(coefficient, se, data=neigh.data[which(neigh.data$boushey==2),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(mo.neigh)

g.neigh <- metagen(coefficient, se, data=neigh.data[which(neigh.data$boushey==3),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(g.neigh)

## Ideological Distance
r.gncp <- metagen(coefficient, se, data=gncp.data[which(gncp.data$boushey==1),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(r.gncp)

mo.gncp <- metagen(coefficient, se, data=gncp.data[which(gncp.data$boushey==2),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(mo.gncp)

##Federal Intervention
r.fed<- metagen(coefficient, se, data=fed.data[which(fed.data$boushey==1),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(r.fed)

mo.fed<- metagen(coefficient, se, data=fed.data[which(fed.data$boushey==2),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(mo.fed)

g.fed<- metagen(coefficient, se, data=fed.data[which(fed.data$boushey==3),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(g.fed)

## Per Capita Income
r.percap<- metagen(coefficient, se, data=percap.data[which(percap.data$boushey==1),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(r.percap)

mo.percap<- metagen(coefficient, se, data=percap.data[which(percap.data$boushey==2),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(mo.percap)

g.percap<- metagen(coefficient, se, data=percap.data[which(percap.data$boushey==3),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(g.percap)

## Legislative Professionalism
r.legprof<- metagen(coefficient, se, data=legprof.data[which(legprof.data$boushey==1),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(r.legprof)

mo.legprof<- metagen(coefficient, se, data=legprof.data[which(legprof.data$boushey==2),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(mo.legprof)

g.legprof<- metagen(coefficient, se, data=legprof.data[which(legprof.data$boushey==3),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(g.legprof)

## Citizen Liberalism
r.cit.lib <- metagen(coefficient, se, data=cit.lib[which(cit.lib$boushey==1),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(r.cit.lib)

mo.cit.lib <- metagen(coefficient, se, data=cit.lib[which(cit.lib$boushey==2),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(mo.cit.lib)

g.cit.lib <- metagen(coefficient, se, data=cit.lib[which(cit.lib$boushey==3),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(g.cit.lib)

## Government Liberalism

r.gov.lib<- metagen(coefficient, se, data=gov.lib[which(gov.lib$boushey==1),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(r.gov.lib)

mo.gov.lib<- metagen(coefficient, se, data=gov.lib[which(gov.lib$boushey==2),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(mo.gov.lib)

g.gov.lib<- metagen(coefficient, se, data=gov.lib[which(gov.lib$boushey==3),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(g.gov.lib)

## Initiative Available
r.init<- metagen(coefficient, se, data=init.data[which(init.data$boushey==1),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(r.init)

mo.init<- metagen(coefficient, se, data=init.data[which(init.data$boushey==2),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(mo.init)

g.init<- metagen(coefficient, se, data=init.data[which(init.data$boushey==3),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(g.init)

## Democratic Control of Government
r.dem<- metagen(coefficient, se, data=dem.data[which(dem.data$boushey==1),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(r.dem)

mo.dem<- metagen(coefficient, se, data=dem.data[which(dem.data$boushey==2),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(mo.dem)

## Republican Control of Government
r.rep<- metagen(coefficient, se, data=rep.data[which(rep.data$boushey==1),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(r.rep)

mo.rep<- metagen(coefficient, se, data=rep.data[which(rep.data$boushey==2),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(mo.rep)


## Divided Government
r.div <- metagen(coefficient, se, data=div.data[which(div.data$boushey==1),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(r.div)

mo.div <- metagen(coefficient, se, data=div.data[which(div.data$boushey==2),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(mo.div)

g.div <- metagen(coefficient, se, data=div.data[which(div.data$boushey==3),], studlab=authors, comb.fixed=FALSE, comb.random=TRUE, prediction=TRUE, sm="OR", method.tau="PM", hakn=TRUE)
summary(g.div)

##################################################################
########################## Figure 4 ##############################
##################################################################

jpeg("figure4.jpg", height=8, width=8, units="in", res=600)
par(mfrow=c(2,2))
funnel(m.neigh, main="(a) Neighbor Adoptions")
funnel(m.gncp, main="(b) Ideological Distance")
funnel(m.init, main="(c) Initiative Available")
funnel(m.rep, main="(d) Republican Control of Government")
dev.off()

funnel(m.fed)
funnel(m.percap)
funnel(m.legprof)
funnel(m.cit.lib)
funnel(m.gov.lib)
funnel(m.init)
funnel(m.dem)
funnel(m.div)

## Reproduction Script for 
## Mallinson, Daniel J. (2020). "Policy Diffusion Results: A New Database for Systematic Review and Meta-Analysis
## State and Local Government Review 52(1): 18-27

## Database Paper Figure 1
tiff("slgr_figure1.tiff", height=6, width=6, units="in", res=600)
funnel(m.legprof, main="Legislative Professionalism")
dev.off()

## Database Table 2
length(unique(meta$model_id[meta$neigh==1]))
length(unique(meta$model_id[meta$geog==1]))
length(unique(meta$model_id[meta$south==1]))
length(unique(meta$model_id[meta$fed==1]))
length(unique(meta$model_id[meta$ideology==1]))
length(unique(meta$model_id[meta$berry_govt==1]))
length(unique(meta$model_id[meta$berry_cit==1]))
length(unique(meta$model_id[meta$ideological_dist==1]))
length(unique(meta$model_id[meta$duration==1]))
length(unique(meta$model_id[meta$slack==1]))
length(unique(meta$model_id[meta$percap==1]))
length(unique(meta$model_id[meta$demog==1]))
length(unique(meta$model_id[meta$urban==1]))
length(unique(meta$model_id[meta$population==1]))
length(unique(meta$model_id[meta$political==1]))
length(unique(meta$model_id[meta$interest_groups==1]))
length(unique(meta$model_id[meta$initiative_avail==1]))
length(unique(meta$model_id[meta$dem_control==1]))
length(unique(meta$model_id[meta$rep_control==1]))
length(unique(meta$model_id[meta$divided==1]))
length(unique(meta$model_id[meta$profession==1]))
length(unique(meta$model_id[meta$context==1]))
length(unique(meta$model_id[meta$other==1]))
length(unique(meta$model_id[meta$intercept==1]))
